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ABSTRACT 






The basic principles of the Galerhin finite element 
method are discussed and applied to two different 
formulations; one using different basis functions and the 
other using the vorticity — divergence form of the shallow 
water ecuations. lach formulation is compared to the 
primitive form of the equations developed by Selley (1976). 
The testing involves a comparison of three finite element 
prediction models using variable size elements. Equilateral 
elements significantly improve the solution and are used in 
most of the comparisons. The formulation using different 
basis functions- produces poorer results than the primitive 
formulation. The vorticity-di vergence formulation produces 
superior results- while eiecuting“faster than the primitive 
model. Fcwever, it does require more storage and the 
relaxation parameters are sensitive to the domain geometry. 
The computer implementation for the vorticity-divergence 
model is discussed and the source listing is included. 
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I. INTROEUCTION 



Shv.man '19?S) claims that progress in numerical modeling 
of the general circulation has heen to some degree dictated 
in the past hy the rate of development in the field of 
computer technology. However, the limited ability to 
oaramet er ize the effects of small-scale processes in terms 
of large scale motions has been an equally important 
limiting factor. Essentially, the major problem of numerical 
modeling of the general circulation is simply that of 
producing a very long range numerical weather forecast. 

Certainly the equations used in the models must be more 
sophisticated to include those physical processes- which are 
u.nlmpor tant for a short range forecast, but may become 
crucial as the length of the forecast is extended. Another 
area where concentrated efforts have improved the forecast 
involves the computational techniques employed to 
approximate and solve the governing equations of the models. 

The motivation behind this thesis is to investigate the 
application of a relatively new computational technique to 
the field of numerical weather prediction. The finite 
element method, long established in engineering, has been 
seriously considered only during the past decade in 
meteorology, -his method has great potential for application 
in atmospheric prediction models. 
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A. BACKGROUND 

The irost coirircn nurrerical integration procedure for 
weather prediction has ^een the finite difference .Tethod in 
which the derivatives in the differential equations of 
notion are replaced by finite difference approximations at a 
discrete set of points in space and time. The resulting set 
of equations, with appropriate restrictions, can then be 
solved by algebraic methods. Until recently, the finite 
difference method has been the workhorse in atmospheric 
prediction models, from their first computer implementation 
to the present . 

’vith the introduction of each new generation of 
computers. the gap between numerical forecasts and 
atmospheric observations has decreased. The rate at which 
this gap decreased has slowed down and appears to be 
leveling off. This would indicate that computer technology 
may not be the primary obstruction to better numerical 
forecasts. In fact, bigger and faster computers alone have 
demonstrated their inability to significantly improve the 
numerical forecast. 

?or example, a major limiting factor of finite 
difference approximations is the truncation error. The 
National Weather Service 7 Layer Primitive Equation model 
( 7L?E Model), operational from 1966 to 1980, had truncation 
errors which increased at a rate proportional to the square 
0 *' the grid spacing. That is, the smaller the grid interval. 
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the smaller the truncation error. To increase its accuracy 
would require increasing the grid iratrix density. This would 
require increased cotrputer storage and corrputat icnal tiire. 
State of the art corrputers are capable of providing these 
additional resources. 

The problem now goes beyond numerical techniques and 
computer technology. Operationally, the National Weather 
service is not capable (due to monetary restrictions) of 
providing a denser concentration of atmospheric 
observations. Therefore, with the present density of initial 
data (observations) and objective analysis techniques 
(getting the data for grid points by Interpcl atlng from 
observed data sources), reducing the grid spacing further on 
the 7LPZ Model does not significantly Increase the accuracy 
c f the solution . 

This additional computer capability can not be utilized 
usins finite difference methods. Therefore, new numerical 
integration techniques must be investigated, such that given 
the same density of observed data, superior sclutlons are 
produced . 

Two alternative techniques, the spectral method and the 
finite element method, have started to gain attention. Toth 
the spectral and finite element methods require more 
computational time per forecast time step than does the 
finite difference method. For example, the finite element 
method requires an equation sclver to invert a larger matrix 
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at each tlrre step for each variable. In this sense, these 
Tethods were held hack hy computer technology, hut recent 
advances in computer technology (i.e. larger and faster 
storage devices, multi-processors, etc) have made these 
alternative numerical techniques competitive. 

;or long range weather predictions, the spectral method 
applied over the glohe or hemisphere is a natural method, 
due to the existence of efficient transforms for the 
nonlinear terms on spherical geometry. It also eliminates 
the truncation error for the horizontal space derivatives 
and the nonlinear instability (aliasing). For these reasons, 
global spectral models have been developed and implemented 
on an operational level, replacing the global finite 
difference models. 

However, because the spectral harmonics are globally 
rather than locally defined, it is thought that for problems 
of more detailed limited area forecasting, the finite 
element method is more suitable. Pioneering work to adapt 
finite element methods to meteorological applications has 
been done by Cullen (1973,1974 and 1979), Staniforth and 
Mitchell ^1977), Hinsman (1975) and lelley (1976). The most 
recent finite element meteorological model at the Naval 
Postgraduate School was written by Kelley (1976) with the 
collaboration of Tr. R.T. Williams. It is this study that 
will serve as a basis for this thesis. The model written by 
Kelley will be altered and used for comparative testing with 
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irproved finite elerrent forrrs implemented by this author. 
Some of the techniques and codes developed by Kelley are 
alsc employed in this thesis. Older (1981) developed a 
technique to smoothly vary the grid geometry in the domain. 
This technique is also implemented both on Kelley's model 
and with the new formulation to give greater versatility 
when testing the model performance. 

3. C3J2CTI7KS 

The objectives for this thesis can be divided into two 
categories: 1) meteorology, 2) computer science. First, the 
meteorological objectives of developing improved finite 
element forms for shallow water equations are as follows: 

1) - Older (1981) after collaboration with Dr. 
M.J.?. Cullen, showed how equilaterally shaped elements 
produced significantly better results than did other 
triangular elements. Kelley (1976) used right triangular 
elements in the implementation of a two-dlmension.al finite 
element model using the primitive form of the shallow water 
equations. A considerable amount of small-scale noise was 
observed in the solution. Hereafter, this model, which was 
developed by Kelley (1976), will be referred to as the 
primitive model. This first objective involves 
re-implementing the primitive model using equilaterally 
shaped elements and comparing the results to those in 
Kelley's thesis. 
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2^ - '•i'illlaT'S and Zieniciewicz (1981) presented new 
finite elerrent techniques for forirulations for the shallow 
water eorations, which use differently shaped functions to 
epprciirate the different dependent variables, which in 
effect stagger the variables. Schoenstadt (1980) 
detronstrated the advantage of spatial staggering of 
dependent variables in finite difference models. The 
application of this technique to finite element models is a 
natural extension, and excellent results were obtained by 
Williams and Zienkiewicz (1981) from application of these 
new formulations on linearized one dimensional cases. The 
objective here is to implement the new forms on the 
primitive model and again do quanltative comparisons of the 
results . 

3) - The major emphasis in this study deals with the 
implement ation and comparison of the vorticity divergence 
form of the shallow water equations, which is described in 
detail in Chapter III. This formulation has the following 
advantages. First, the geostrophlc adjustment process is 
treated better than in the primitive form of the equations. 
Secondly, the velocity and height fields are evaluated at 
the same grid point, where the best primitive form requires 
staggering these dependent variables. And thirdly, a larger 
time step is allowed due to the semi implicit form of the 
equation. Again comparisons between the results from the 
vorticity divergence and primitive model are presented. 
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The computer science aspect of this thesis was primarily 
devoted to the implerrentat ion cf the different models and 
the style and architecture of the program. Finite element 
methods require more computational time than do finite 
difference methods, not only in the solution of the 
equations, hut also in the amount of computation required to 
evaluate each term in the equations. 

The implementations of these two dimensional models, 
although complex when viewed from the surface, have a lot of 
generality and redundancy in the operations required. 
Versatile modules can he written to ease the implementation 
and facilitate changes. The objective here is to efficiently 
implement these new forms and demonstrate the utility of 
these versatile modules for future implementations. 

C. TEISIS STRUCTURE 

This thesis presents the results obtained from tests of 
the various finite element formulations. The resul-ts are 
compared to those from the primitive model written by lelley 
(1976). Accompanying the results is a detailed discussion of 
the reformulation and implementation process. 

Chapter II of the thesis presents a tutorial of the 
finite element method and the area coordinates system used 
in the evaluation of the element integration. The Galerkin 
finite element method used in all the models is developed 
and applied to the advection equation in one dimension. 
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Cicpter III presents the detailed description of the 
VC r i ci ty-di vergence shallow water .T’odel. Here the equations 
are shewn and written using the ^alerkin .Tethod. A 
discussion of the corputat ional technique used is presented 
along with the model's physical parameters. 

Chapter IV presents a descriptive overview of the 
computer implementation. The chapter includes a list of 
options available for testing, a brief description of the 
matrix compaction technio.ue and the formulations of the 
versatile modules used to implement the complex equations. 

Chapters V through VII discuss the results obtained from 
the different experiments. Chapter V briefly describes the 
primitive model used for all comparisons and the results 
*'rom changing the element shape to equilateral triangles. 
Chapter VI reformulates the primitive model so that the 
geopotential is staggered with respect to the velocity 
variable. Jor simplicity, the continuity equation is also 
linearized. Chapter VII compares the results from the 
vort i ci ty-di vergence model developed in Chapter III to those 
from the primitive model. 

The last chapter summarizes the results from all the 
experiments and identifies what areas require follow on 
work. The source code for the vorticity-divergence model is 
presented in Appendix A. 
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II. FINITE SLZf^ENTS 



As is often the case with an original de velopr'ent , it is 
rather difficult to quote an exact date on which the finite 
elerrent trethod was invented, but the roots of the method can 
be traced back to these separate groups: applied 
rratherraticians , physicists and engineers. Since the early 
d eveloprrents of the finite elerrent rrethod, a large amount of 



research has been devoted 


to the 


technique 


However, 


the 


finite element method 


obtained 


its 


rea 1 


impetus by 


the 


independent developments 


carried 


out 


by 


engineers. 


Its 



essential simplicity in both physical interpretation and 
mathematical form has undoubtedly been as much- behind its 
popularity as is the digital computer which today permits a 
realistic solution of even the most complex situations. 

The name " finite element *’ was coined in a paper by 
?..V. Clough, in which the technique was presented for plane 
stress analysis, as discussed in Bathe (ISVe). While finite 
element methods have made a deep impact via the field of 
solid mechanics, where it can be said that today they 
represent the generally accepted method of discretizing 
continuum problems for computer-based solution, the same 
appears not to be true in fluid mechanics or atmospheric 
pred iction . 
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.\urrerous finite element formulations are currently 
availatle. Strang (1973), Ncrrie (1373) and Zienkiewicz 
'1971' present detailed theoretical discussions of each, -he 
C-alerkin method, the most popular finite element method, is 
described in detail below and used in the equation 
formulation later. 

A. FINITI ZLZr^INT CON’CZPT 

-he problem of solving partial differential equations 
can be specified in one of two ways. In the first, finite 
difference methods specify the dependent variables at 
certain grid points in space and time, and the derivatives 
are evaluated using Taylor series approximations. Secondly, 
the calculus of variation requires the minimization of a 
functional over a domain, where a functional is defined as a 
variational integral over the domain. 

The calculus of variation approach creates a purely 
physical model where the functional equivalent to the known 
differential equations are known. Its major disadvantage is 
that it limits the method only to those problems for which 
functionals exist. Finite element methods, an extention of 
this method, derive mathematical approxima t ion s directly 
from the differential equations governing the problem. The 
advantage here is that it extends the method to a range of 
problems for which a functional may not exist, or has not 



beer, discovered. 
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The finite elem 


ent 


method 


di vides 


the 


domain into 


subdomains or finite 


elements (u 


sually of 


the 


same form). 


Nodes are located 


along 


the 


boundary 


of 


the elements. 



usually at the element vertices and at strategic positions 
fmidside, centroid, etc.) in the interior and on the sides 
of faces of eleir.ents. 

CofTmonly used elements are triangular, polygonal or 
polyhedral in form for two-dimensional problems. The choice 
of elements depends on the type of problem, the number of 
elements desired, the accuracy required and the available 
computing time. To begin with, the element must be able to 
represent derivatives of up to the order required in the 
solution procedure, and to guarantee continuous first 
derivatives across the element boundaries to avoid 
singularities. Triangular elements are employed in this 
thesis because they can be used effectively to represent 
irregular boundaries, and/or geometry, and also to 
concentrate coordinate functions in those regions of the 
domain where rapidly varying solutions are anticipated. 

Consider the problem of solving approximately the 
differential equation 

L(u) = f(i) II-l 

where L is a differential operator, u the dependent 
variable, and fix) is a specified forcing function. Suppose 
that II-l is to be solved in the domain a ^ x - b and that 
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appropriate boundary conditions are provided. The residual R 
is forred fror II-l as follows: 

L(u) - f{x) = R I I -2 

The critical step is to select a trial farrily cf 
approximate solutions (the members of a trial family are 
often called basis functions). The basis function is 
prescribed 'functionally) over the domain in a piecewise 
fashion, element by element, and are generally a combination 
of low order polynominal s . A one dimensional example is 
shown in Figure 1, wherein the domain (i axis) is divided 
into six elements (line segments) A through F. The basis 
functions are linear and one is shown for node 4 only in 
Figure 1. The function has a value of unity over node 4, and 



decreases linearly 


to is zero at nodes 3 and 5 


and 


zero 


els ewhere. 












Consider 


a 


series 


of linearly independent 


basi s 


functions V. 

J 


(i) , 


as in 


Figure 1 . Now u (x ) 


.can 


be 


approximated 


wi th 


a finite 


series as follows: 








u(x 


) = I 
1 ^ 


V. (x) = V. 

J J J 




II -3 



where i.is the coefficient of the jth basis function and has 
a value equal to u at node j. 

Substituting this approximate solution II-3 wherever u 
appears in the differential equation II-l 

L(i^.V.) - f(x) = R 1 1 -4 

J «J 



22 




23 



The best solution will he one which in sooie sense 
reduces the residual R to a minimum value at all points in 
the domain. 3y definition, the residual obtained using the 
exact differential equation is identically zero everywhere. 
The residual R, formed in eo.uation II-4, is minimizied when 
multiplied with a weighting function, integrated over the 
domain and set equal to zero. This process is known as the 
weighted residual method 



where '* is the weighting function and is referred to as the 
test function in the following development. The weighted 
residual method minimizes the errors of the residuals, such 
that the summation of all the positive and negative errors 
add to zero . 

The Galerkin method, the most popular fi ni te .element 
method, is more general in application and is a special case 
of the method of weighted residuals, as discussed by Finder 
and Gray ^1977). The requirement imposed on the weighted 
residual method forming the Galerkin method is: 

* the test (weighting) function be equal to the 



b 




II-5 



a 



basis (trial) function W = V . This process 
leads in general to the best approximation of 



the solution. 
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The final ^alerkin form is obtained by substituting II-4 
into II-5, yielding 



b 

J 

a 



b 



{ X )dz = 2 

a 



1 1 -6 



If this procedure is repeated for N points in the dorrain 
a system with N equations and iV unknowns will be generated. 



3. galepsin application 

The following ezample taken in part from Haltiner and 
Williatrs (1S80) applies the Galerkin method to the advection 
equation with linear elements 

5u Bu 

— + c — =0 II-7 

3t dx 

This equation is dependent in both time and space. The 
treatment of time variation is important. for most 

meteorological prediction problems. The Galerkin method is 
not applied to the time dependence because it is more 
convenient to use finite differences in time, as is done with 
this example later. The same treatment is applied tc the 
prognostic equations later, where two finite differencing 
methods are employed to do the time integration. 

The Galerkin procedure represents the dependent variable 
u(x.t) with a sum of functions that have the prescribed 
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spatial structure as in Figure 1. Approrirrate u(i,t) with 
the finite series as follows 



u ( X , t ) 






1 ^At)V.(x) 
J=l ^ 




II-8 



where the coefficient ^ . (t), a function of tixe, is the 

tJ 

scalar value of u at node j. The basis functions, V.(x), are 
functions of space only and j equals 1 to 7 for the exatrple 
in Figure 1. The repeated subscript in this forx implies a 
surr ever the repeated subscript. 

The Galerkin equation for the advection equation II-7 is 
obtained by setting L = c(d()/dx) and substituting in the 
approximate solution II-8 wherever u is found. 



b b 

K r W r dV. 

1 V.?. dx + c I 0. — = 0 II-9 

j-i at J J ^ j=i J ax ^ 

a a 



where i = 1 to N, the test function and ^he basis 

function. The domain of integration is given by a - x - b, 
and the integration is done in a piecewise fashion, element 
by element . 

In this one-dimensional case, an equation like II-9 is 
written for each node, i. Considering node 4, what are the 
possible non-zero contributions from equation II-9? Figure 
2 illustrates the basis and the test function interaction 
during the piecewise integration process. From the 
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1 



Basis Function 



Test Function 
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Figure 2 . Basis aind test function interaction 
during the piecewise integration process. 
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'iefir.ition of the basis and the test function, locally- 
defined as unity at node j and linearly decreasing to zero 
-t J ± 1 and zero elsewhere, the only non— zero 
contributions are made when J = 3 over element C, j = 4 over 
elements " and I, and J = 5 over element D. 

The evaluation of II-9 for i = m, which is given in 
Haltiner and Williams fi981), leads to the equation; 



i ±. 
e dt 
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0 11-10 



The boundary points, which in this example are nodes 1 
and 7, are evaluated in the same way as the interior nodes, 
with the exception that cyclic conditions are imposed. 

The time discretization of'-’I I-'i0“ls d-one using- a " fl’nite 
difference scheme. Applying leapfrog time differencing gives 
the following equation 
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The resultant equation set, in matrix form, contains an 
N’xN matrix where N is the number of nodes. 

The transition from one-dimension to two is 
mathematically Identical. The domain is now subdivided into 
finite areas, which are triangles in this implementation and 
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the basis functions are linear. However, now they are 
pyrarrid shaped with value unity at the center and decrease 
to zero at the surrounding nodes, and are zero elsewhere. 
Figure 3 shows this basis function for node 28 outlined in 
heavy black. The value at any node again can be approximated 
by II-3, where j ranges over all nodes connected to node i 
including i itself. The connectivity for node i = 28 in 
Figure 3 is j = 15,16,27,28,29,39 and 40. 

The integration is still over the entire domain. Vith 
both the basis and the test function zero over the domain, 
except locally over each element, the global integration can 
be performed by integrating locally over each element. By 
definition, this integration can be expressed as an inner 
procuct of both functions (i.e. basis, test) as follows: 



Using this definition and the repeated subscript 
notation equation II-S becomes 



where the dot implies differentiation with respect to time, 
and the second subscript implies differentiation with 
respect to the second subscript. The local Integration may 
be calculated directly from exact expressions derived from 
area coordinates described in detail in the next subsection. 
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Figure 3* Basis function for node 28. The 
shaded area is the complete tasis function 
and the V. , where j = 15, 16,27,28,29,39,^0 
are jth node basis fxinctions for node 28. The 
dashed line at node 28 has length unity. 



30 



In surrrrary, the Galerlcln procedure Involves subdividing 
the domain into finite elements, approximating the dependent 
variables by a linear combination of low order polynomials 
and substituting them into the equations. The equation is 
multiplied by a test function, integrated over the entire 
domain and finally the resulting system of equations is 
solved . 

C. AREA COORri NATES 

■»’hile the Cartesian coordinate system is the natural 
choice of coordinates for m.ost two dimensional problems, it 
is not convenient when working with triangularly shaped 
elements. It is therefore necessary to define a special set 
of normalized coordinates for a triangle. Area, or natural 
coordinates as they are commonly called, reduce the 
formidable task of Integrating products between the basis 
and test functions and their derivatives over a triangular 
element and result in easily -eomputable and exact 
expressions. 

The following development is taken in part from the 
formulation by Zienkiewlcs (1971). Consider the triangular 
element Illustrated in Figure 4. There is a one-to-one 
correspondence between the Cartesian coordinates (X,Y) and 
the area coordinates ) for the element. Let A 
denote the area of the triangle and Aj^ , ^and A^ the areas of 
the subtriangles in Figure 4 such that A = Aj^+A^-^A^. 
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( 0 , 0 . 1 ) 




Fig 4. Cartesicin vs. area cooirdinates 
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The relationship between a point ?(X,Y) In Cartesian 
coordinates and L^, ) in area coordinates can be 

seen tv the following transforrratlons 

X = Lj_X + L2X ^ L^X 

Y = L^Y + L 2 Y + L 3 Y 11-14 




and 

= (2A bj_X + aj_Y) /2A 

L 2 = (2A + b 2 X a 2 Y) /2A 11-15 

= ^2A + b^X ^ a^Y) /2A 

where 2A Is twice the area of the triangle and the a's and 
b 's are defined as In Figure 5. 

It is worth noting that every tuple (L , ^2 ’ ^3 ^ 

corresponds to a unique pair (X,Y) of Cartesian coordinates. 

In Figure 4, Lj^= 1 at vertex 1 and 0 at vertices 2 and 3. A 
linear relation exists between the area and Cartesian 
coordinates which Implies that values for L^^vary linearly 
over the triangle with a value one at vertex 1 and a value 
of zero at vertices 2 and 3; and similarly for L 2 and . 
This demonstrates how each component in the tuple L 2 1 L^) 

behaves over the triangle as do the linear basis and test 

functions over the element, as was seen in Figure 4. Clearly 

= Vi 11-16 
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where is a linear function of the Cartesian coordinates 
( 1 . e . basis, test ) . 

Zienkiewicz (1971) shows that it is possible to 
integrate any polynorrial in area coordinates using the 
sirple relationship 
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where rr , n and p are positive integers and A is the 
elen'entary area. For an example of this integration 
technique using inner product notation, equation 11-12 is 
evaluated as follows 
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The differential operations in area coordinates follow 
directly from the differentiation of (II-15) where 



B 3 b 
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bx 



i=l 2A bL, 



and 



b 3 a^^ b 

by i=l 2A bL, 



11-20 
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flS explained earlier (see Equation 11-16), Vj_is a linear 
function ^i.e. basis, test) which equals a component of 
the area coordinate tuple. Therefore 
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Consequently 







forj=l is 
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As an example, 
vertices j = 2, i = 1 



consider the inner product <V. ,V.> at 

X 

This integration is evaluated as 
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dxdy 
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b_ 1! 0! 0! b, 

= -X 2A = 

2A (1+0+0+2)! 6 

Therefore any inner product in the formulation can be 
readily evaluated using area coordinates. Another benefit of 
using this coordinate system is that all of the inner 
products are functions of space only and need be computed 
only once. 
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III. SHALLOW WATIR MCEZL 



The governing equations for this irodel are derived hy 
rrcki.ng several sir’plifyiag assujptions on the primitive 
ecuations of motion, vhioh then give the barctropic shallow 
water equations. However, as mentioned previously, the 
shallow water equations describe many significant features 
of the large-scale motion of the atmosphere, and therefore 
have been used in numerous experirrents over the years. 

The vcrtici ty-divergence form of the equations has 
several advantages. Williams (1981) has shown that the 
geostrophic adjustment process is treated much better with 
the vcrticity divergence formulation than with a direct 
treatment of the primitive form of the shallow water 
equations, such as was used by Kelley (1976). This 
formulation also allows the- velocity components and the 
height to be carried at the same nodal points, whereas the 
best scheme for the primitive form of the equations requires 
staggering of the fields, as seen in Schoenstadt (1980). The 
vo^ticity divergence form of the equations is also 
convenient for the application of semi-implicit 
differencing, which saves considerable computer time. 
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GOVERNING ICUATIONS 
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Equation 'III-l) Is the continuity equation and the III-2 
and III-3 are the momentum equations, respectively. The 
variables are defined as follows: 

T. ,j - the spatial coordinates of the domain 
u,v - components of the wind vector 

1 geopotential = (gravity x free surface height) 

2 2 

- mean §eopotential = 49,000 meters /seconds 
t time 

2 - kinetic energy 

C - absolute vorticity = ("S + fj 
S relative vorticity 

f, - ccriolis force (mid-channel f-plane) 
r divergence 

The shallow water equations can be written in 
vorticity divergence form as follows: 
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bt 



I II -4 



X a b 

+ C(j) = - — _ — (^v) 
bx By 

BS B B 

— (uQ) (vC) III-5 

Bt Bi by 

2 . ^ ^ 2 

— + (vC) (uQ) - III-6 

Bt Bx By 

where III-4 is the sere continuity ecuation as III-l, III-5 
is the vorticity equation and III-6 is the divert^ence 
equa t ion . 

Because of the vorticity divergence forrr of the 
equations, it becomes necessary to solve the time dependent 
variables "S and I in terms of 4/, the stream function 
(rotational part of the wind), and'X, the velocity potential 
fdivergent part of the wind). The initial fields for the 
model will be in terms of 4(, tC and 0. 

The following diagnostic relationships are defined and 
used later in the solution of the equation set. 



u = - 'Jy* *x- 


III-7 




III-8 


where the subscript implies differentiation 

2 2 
u V 


K = kinetic energy. 


III-9 


uC = u ('S + fj , 


III-10 


vC = v(5 + fj , 


III-ll 


« = 0u , 


III 12 


p = ;z)v, 


III-13 
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I 



'S 



III-14 



C = 'X. III-15 

3. ZCUATION rOP.f'ULATION 

The G-dlerkin method described in Chapter II is now 
aoplied to eonations III-4 through HI-15. For ease of 
comprehension, the shorthand inner product notation as in 
II-12 will be used to simplify the equations. The detailed 
Galerkir. formulation will be shown for equation III-7, the u 
component of motion. The method follows directly from the 
example in Chapter II of this thesis, which in turn follows 
in part from Selley (1976) and Ealtiner and Williams (1981). 

Consider equation III-7 and assume that each variable u, 
^ and C( is aoproiimated by 



u 


“ "j 


’'j • 




'i' 




''j • 


iii-ie 


1C 


II 


''j • 





where the repeated subscripts indicate summation over the 
range of the subscript. Substituting these approximate 
solutions into III-7 yields 

u. V . = ~^.V.) + — {‘5C. '/.) 1 1 1-17 

Since only the basis function 7. is a function of space, 

c 

III-17 ray be further simplified "b'/ factoring out the time 



deoendent coefficients. 



The next step requires multiplying by a test function 
as liscussed in Chapter II, and integrating over the area 
dome in 



UjJj V.V^dA = - dA 

A 



ff ^V, 



+ X. 



J y 



— ‘'y. dA 
bx 
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The final form in inner product notation is 






- Vjy.V,> * <5(jVj^.V,> III-19 



where the double subscript implies differentiating with 
respect to the second subscript. 

The three prognostic equations (III-4, III-5 and III-6) 
are similarly advanced using the Galerkin technique to 
become, respectively: 
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<■5.7 ..V,> = -<(u0) .7. ,7. > - <(vC) .7. ,7, >111-21 
Jji Jjxi jjyi 



+ <K . v^7 .,7. > 1 1 1-22 

J J 1 

2 

•rfhere v is the Laplacian operator and the dot implies 
differentiation with respect to the time dependence in 
III-4, III-5 and Ill-e. 

Similarly, Galerkin equations are formulated for 
ecuations III-7 through III-15. 
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C. Zli'.I riSCRSTIZATION 

The eauation set III-20, III-21 and III-22 is arranged 
so that all the terrrs on the left hand side can be treated 
ijplicitlv, and all the terms on the right hand side can be 
treated explicitly. The explicit time integration will be 
done by the leapfrog difference method. To start the time 
integration, two forward half steps are taken, after which 
the full leapfrog scheme is used for the remainder of the 
forecast period. 

The vorticity equation III-21 is solved Independently 
from III-20 or III-22. However, III-20 and III-22 
^continuity and divergence equations, respectively) are 
coupled. To explicitly solve either, decoupling of the 
equations is necessary. In this thesis this is done through 
algebraic substitution of III-22 (solved for D(n+1)) into 

the time integration is performed on HI-20, 
solved for D(n + 1) using the jj (n+1) value, 
prediction equations are 

C<Vj,V^>] 

- - [5EP.YJ 

- . <Vjy.Vjy>] 

- 3[(^u)5<Vj^.Ti > * (0v)5<Vjy.Vj >] 

III-23 



HI-20. Cnee 
HI-22 can be 
The final 

^ j jx * ^ ix^ 



1 



whs-e A = 4/^2At) . B = A/^ , C = B/(2dt) and [BCRY] is the 
ssostrophic boundary contribution, see Section 



■Sj <Vj.vi> 












in-24 



tr-'b-’i 



r"J^<Vj,Vi> - (At/2) j,v^> 



III-25 



After these three elliptic equations are solved, the 
history of the variables III-7 through HI-15 is updated. 

A large tirre step can be applied to this fortr of the 
shallow water equations due to the semi-i.^plicit nature of 
the equations. "Ihis is very important since finite element 
methods generally require more computer time per time step. 
The VC r ti ci ty-d i vergence formulation acts as a filter, which 
slows down the high frequency waves in the solution. The 
two-dimensional advective stability criterion for a linear 
element, derived by Cullen (1973), was used to determine the 
correct time step, 

a X 

at = III-26 

Id 1^6 

whereat is the time step in seconds, ax the shortest grid 
spacing in meters and c the fastest phase velocity. 
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D. COr^PUTATlONAL TSCHNIQUZS 

The final prognostic equation set requires the solution 
of a I-elrr.holtz equation for ^ and Poisson equations for 'i' 
aad 1. The frost cofriron method of solution used hy 
meteorologists has been the successive over relaxation 
method fSOR) in which an initial guess of the solution is 
made and then progressively improved until an acceptable 

level of accuracy is reached. SOR is employed in the 
solution of the equations, where III-23 can be represented by 

= {b} III-27 

and III 24, III 25 by 

= {b} III-28 

2 

where v the Laplaclan operator, [M] = <V.,V. > matrix, {x} 

o ^ 

- the dependent variable in vector notation, C - constant as 
in III-23 and {b} the right hand side of the equation or the 
forcing function. 

The mass matrix [M] , dimensioned (nxn), is a matrix of 
coefficients whose rows are the equations of the system to 
be solved. There exists a one to one correspondence between 
the rows of the mass matrix and the nodes of the domain. 
Zach equation has a term (column) for each node, where a 
non-zero term represents connectivity. Nodes are connected 
if they are both vertices of the same element. Obviously [M] 
is a sparse matrix containing the inner products for the 
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left hand side. Chapter 4 of this thesis will describe the 
ratrix coirpaction procedure. 

The forcing function {h}, dirrensioned (nxl), involves 
only variables at the current time step and is easily 
computed using four very versatile subroutines described in 
detail in the next chapter. 

The initial guess to start SCR is the previous time step 
solution. An average of 30 passes per equation are needed 
for each time step. The solution is considered to have 
converged to its final value when the residual for each node 
has been reduced to some acceptably small value. 

The diagnostic equations III-7 through III-15 must also 
be solved every time step. Eowever, the same technique is 
not used for these equations. Dr. Y.J.P. Cullen suggested an 
under relaxa ti on scheme for which three passes over the 
dcrain should produce a solution of acceptable accuracy, 
since the coefficient matrix is so strongly diagonally 
dominant, -ass lumping of the coefficient, matrix is used for 
the first guess. This technique requires replacing the trass 
matrix Cm] by the identity matrix [1]. A first guess of this 
type is able to describe most of the large scale features, 
which in turn reduces the number of iterative passes over 
the field. Successive passes converge to solutions which 
describe smaller scale motion, approximately to the same 
order of magnitude as introduced by computational error, so 
that further iterations are not needed. 
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I. GRID GIOMITHT 

The dotralr. of this rrodel is a cylindrical channel, with 
total length of 4045 Km and width of 3503 Km. The channel 
sirulates a belt around the earth and it proves to be an 
excellent test bed for comparing with the finite element 
formulations used by Kelley (1976) and Older (1981). 

The domain is subdivided into’ equilateral triangles as 
shown in Figure 6. I^ost of the test runs for this thesis use 
a 12x12 mesh which has 156 nodes and 288 elements. This 
Implementation is not restricted to one grid pattern. The 
technique developed by Older (1981) to vary the nodal 
geometry smoothly to achieve areas of denser and coarser 
resolution is also implemented, as in a third grid pattern 
that varies the nodal geometry abruptly. A short discussion 
of these nodal geometries with accompaning illustrations of 
each is presented in Chapter VII, where the different test 
cases are described. 

Cyclic continuity is assured in the x direction by 
wrapping the domain around the earth to form a cylindrical 
domain. This has the advantage of eliminating the east-west 
boundaries and it simulates the flow around the earth. The 
only boundaries cn this domain are the north-south walls and 
their treatment will be discussed shortly. 
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Figure 6 Domain divided into equilateral triangles. 
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INITIAL CONDITIONS 



As rreationed previously, the ref orirula t Ion of the 

scverning equations into the vorticity-diver^’ence 
shallow water equation set requires solving the tirre 

dependent variables in terms of the stream function and 
velocity potential. The continuity equation is not altered, 
so that its solution is expressed in terms of 0. 

For the basic testing of the model's performance, simple 
analytic sinusoidal initial conditions are used to insure 
the most accurate analysis possible and to simplify the 
compari so ns • 

The sinusoidal initial fields are graphically shown in 
Figure 7 as 3-dimensional surfaces. The geopotential field 0 
consist of a half sice wave in the y direction and a single 

cosine wave in the x direction. The stream function 

calculated by dividing the geopotential field by the 
Coriolis force, has the same physical structure as 0. The 
velocity potential I has a single .sine wave > in- the x 
direction and a half sine wave in the y direction. 

These initial conditions are computed as follows 



0 = fo AsinoCj^cosoc 2 - ^ 

= ^/fo III-29 

1 = Csinocj_sincv 2 quasi-geos t rophic divergence 



where 



A - arbitrary amplitude 

fj - Coriolis value for mid-channel latitude 
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a) 0 and Y initial fields. 




b) X initial field. 



Figure . 3“Gioensional view of the inital fields. 
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U mean flow 

- T.il-lati tude value of y 
f - mean free ^eopotential height 
= 49,200 ,T^/s ^ 

y.2 Zttxt/L 
r - wave nurrter 
V channel width 
L - channel length 
C - -(foUoc2BA)/(f? + $3) 

3 - c<l + o'l 
3. 3CUNEARY CONDITIONS 

Boundary conditions are only required on the north and 
south walls of the grid domain. Due to cyclic continuity, 
the domain is wrapped around creating a cylinder eliminating 
the east and west boundaries. However, careful attention to 
detail is needed during the implementation to assure this 
continuity. Separate boundary conditions are applied to each 
of the predictor equations III-23, III-24 and HI-25. These 
conditions are computed for the wall nodes only and are 
applied during each pass through the relaxation scheme. 

The vorticity eouation III-24, the most sensitive of the 
predictor equations to solve, requires ^ on the north-south 
boundaries to remain constant for the entire forecast 
period. Since this equation is solved in terms of 'If, the 
initial r. orth-south Sf values are saved and assigned to the 
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boundary points after each pass through the relaxation 
subroutine. 

Che oroper boundary condition for the divergence 
equation III-25 would be dX/dn = C. However, for the purpose 
of this study, there is more interest in the sinusoidal 
variation in the y direction and not in the region of the 
walls. Therefore I = 0 is appropriate. 

The continuity equation III-23, the most complex 
predictor equation, requires that there be no mass flux 
through the north-south walls. The geostrophic boundary 
condition 

— = -uf„ III-30 

dy 

is anplied to the north south boundary nodes for the terms 
[3E?T] in equation III-23. Integrating the inner product 
<i^.. V.,V. > by parts produces the boundary terms 

jf v^(J 5.V.)7 dxdy = f/ ^ ^ . V . ) V . dzdy 

yX J J 1 yj J J 

= l( [y v(^jVj ) ) - 7(5^ jV j)»v7 dxdy 
yx 

= fBERI] - ^ III-31 

where n is a unit vector normal to the domain and dr is the 
differential distance along the path of integration on the 
perimeter cf the domain. 
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The geostrophic boundary condition III-32! is substituted 
into the contour integral in equation III-31 and put into 



advective equation in Chapter II. The resulting term is 
derived as follows 



Zquation III-32 appears twice in the continuity equation 
III 23, for time levels (n-^l) and (n-1). All values of u are 
known •'or time ^n-1), since they are saved from the previous 
calculations. However, u(n*l) has not been computed. To 
solve for u(n+l), both 'l'(n + l) and T(n+1) are needed. 4f(n+l) 
is solved first from the vorticity equation, ^(n+l) needs 
as part of its solution and 0(n + l) needs u(n + l) in 
its solution. To avoid this problem, it is assumed that 
7fn+l) has a negligible con tri bution to the solution of 
u(n+l) and only'9f(n+l) is used. 




the same way as in the one-d im.ensional 




(u . , + 2u. * u. , III-32 

3 J+1 J J-1 i 
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IV . CCrPUTZE ir^PIIf^ZNTATION 

-he f or.Tulation and general iheory of the finite ele.Tent 
rethod was presented in the previous chapters. The objective 
ir this chapter is to discuss sore important computational 
aspects pertaining to the in-pleirenta t ion of the finite 
element prediction system. 

The main advantage that the finite element method has 
over other prediction t echniques is its generality. 
Conceptually , it seems possible by using many elements, to 
approximate virtually any surface with complex boundaries 
and initial conditions to such a degree that an accurate 
solution can be obtained. In practice, however, obvious 
engineering limitations arise, a most important one being 
the cost of the computation. As the number of elements is 
increased, a larger amount of computer time is required for 
a forecast. Furthermore, the limitations of the program and 
the computer may prevent the use of a large number of 
elements. These limitations may be due to the computer speed 
and storage availability, or round-off errors propagated in 
the computations because of finite precision arithmetic. 
Also, the malfunction of a hardware component, if the 
prediction is carried out using many computer hours to 
execute, can be a serious problem. It is therefore desirable 
to use efficient finite element programs. 



The effectiveness of a program depends essentially or. 
the folloving factors. Firstly, the use of efficient finite 
elef^e.nt tech-iiques is important. Secondly, efficient 
prc5:ramrri ng methods and sophisticated use of the available 
computer hardware and software are important. The third very 
important aspect in the development of a finite element 
nrosrem is the use of appropriate num.erical techniques. 

The vortici ty-divergence model described in the previous 
chapter is implemented on the 3333 computer located at 
the N'aval Postgraduate School. Some notable features of its 
architecture are the three trillion bytes of virtual mass 
storage, of which eight mega bytes are available to each 
user, and the 57 nanosecond machine cycle time. The model is 
executed mostly using a 12x12 element domain requiring 430h 
bytes of storage and 33 seconds of CPU time to execute. 
Inceedir.g execution time and/or available storage is not a 
problem, in fact the system allowed a lot of flexibility 
during the impl-ementat ion. phase of the model. 

The source code is written using FORTRAN IV and compiled 
on an optimizing Fortran H compiler. Appendix A contains the 
source cede listing, which is divided into five subdivisions 
delineating the logical structure of the program. 

A. PROGRAM ARCEITFCTURZ 

Program features incorporated in the model are: 

1) 'Modularity. With only a few exceptions, each 
module is limited to one page of FORTRAN code. This makes it 
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eesier to comprehend the program. Each module performs only 
one task. Ecr example, subroutine CONTEO computes the value 
of the forcing ter^s for the continuity equation. Likewise 
there is also one module for the divergence and vcrticity 
equations. To implement a new set of equations, only these 
rcdiles would have to be altered. 

2} Easily controllable switches. Switches may be set 
to either print, plot or tabulate harmonic analysis data for 
most of the available fields. The ability to display 
intermediate results allows each portion of the algorithm to 
be monitored for computational adjustments. This also makes 
it easier for unfamiliar users to become acquinated with the 
computational model. 

3^ Eorcing term subroutines. In previous 
implementations, each forcing term was calculated by a 
special subroutine. In this implementation, the calculations 
are accomplished by general purpose routines which simplify 
the implementa-tion. of the compjl.ex prognostic eq-uations-. Th±:S" 
allows implementa tion of different equation sets (i.e. 
^aroclinic ^’odel) over the same domain with minimal effort. 

4) Tccumentation . Each variable is defined by a 
short phrase (Appendix A, A.). The function of each module 
is described in an introductory paragraph. Shared data is 
placed in named common blocks and identified with each 
subroutine which uses them. A subroutine index is given. 



1 • ^aln Prograir 

The main program is short, calling only six modules 
which reflect the basic sequential flow of the model. It 
starts with initialization of all model parameters (i.e. 
model options, domain, finite element arrays, inner 
products). It then initializes the input fields (i.e. 
gecpctential heights, stream function and velocity 
potential) and is followed by initialization of all 
remaining dependent variables. M this time the model is 
totally initialized and time integration begins. As 
mentioned previously two techniques are employed for time 
integration, each having its own module. Upon completion of 
the last forecast, the program terminates. 

Arrays are the only data structures used and are 
srcuped usins IS different common cloclcs. Several arrays are 
used as static link lists, as described in detail later, 
which simplified the algorithms. The common block format has 
the advantage of reducing the overall execut ion- time of the 
program, ^ost of the arguments passed during a call to a 
subroutine are contained in conmon. This requires less time 
to execute since no parameter passing is required for the 
arguments. Another benefit of this format is that the code 
becomes less cumbersome and more readable. Each variable and 
array is defined in the first subsection of Appendix A along 
with a page index for the subroutines. 



c: c; 



Initialisation Phase 



Appendix A, Section C contains all the subroutines 
used durin,£ the initialization pnase of the pro^rarr. From 
the user point of view, the most important suoroutine is 
I.viFSr, the first subroutine called, which contains all the 
global variables that control the different options 
available per run. This is the only subroutine that is 
changed to run the different experiments, assuming that no 
new computational technique is introduced. The selection of 
options are: 

1) - channel location - the channel may be 

shifted north or south by presetting the 
north/south latitude limits in INITG3. 

2) - variable geometry - the node positions may 

be grouped for more dense node patterns to 
yield higher resolution. Two variables HI 
and H2 set the ratio used to vary the 
spacing -along the x .and. y--iaxes, 
respectively. 

3) - initial field wave length and amplitude can 

be altered to produce various effects. 

4) change the initial mean flow. 

5) - diffusion can be entered for any of the 

three prognostic equations. 

e) - maximum length of forecast period may be 
changed and a print, plot or harmonic 



analysis of any dependent variable may be 
requested for any time interval. 

Once the eiperiment is determined, the options 
listed above are set. The program is ready to be executed. 

The largest part of the initialization phase 
consists of establishing the domain and producing all the 
finite element computational vectors that remain constant 
throughout the experiment. 

The first several steps in setting up the domain are 
concerned with indexing. Subroutine COF.RES is called first, 
where all the nodes (grid points) and elements (triangular 
areas) are numbered consecutively starting at the southwest 
corner of the domain and moving eastward across each row or 
latitude. ■ Tor ea'ch elnement, a'^recor d" xf 'all’ 'o f its nodes 
(vertices^ are stored in array (y,3), where M is the 
total number of elements. To facilitate the inner product 
evaluation later, a local numbering system is required for 
each element. That is, fcr each element, its nodes are 
stored counterclockwise in a positive sense. The first node 
however, is arbitrary. 

With the domain divided and numbered, a connectivity 
list (the correspondence between each node and the neighbor 
nodes) is constructed for each node by subroutine CORRZL. 
Each node is adjacent to four or six other nodes depending 
on whether it is a boundary or interior node, respectively. 
These adjacent nodes, plus itself, make up the connectivity 
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list for one node. The connectivity lists are then 
concatenated sequentially starting with the first nodes 
''onnectivity list into the vector NAf^I (NN), where NN is the 
sur of the nodes in each connectivity list. (i.e. for a 
12 t: 12 dorain with 156 nodes, and equilateral eleirents NN = 
lf’44). :5or the first tirre during the initialization phase, 
special attention is given to cyclic continuity. As 
discussed earlier, cyclic continuity is the joining of the 
east and west boundaries to create a cylindrical channel. 
The connectivity list for these east/west boundary nodes 
just be cojplete to insure proper continuity for the 
calculations later. 

The connectivity vector NAf^E is frequently used 
during most computations. Two utility vectors ISTART 
(containing the starting location in NAf'E for a particular 
node^ and NUM (containing the number of nodes in its 
connectivity list) are used to locate and index through the 
vector NAME, as will be seen shortly. This same' technique is 
used to index through the coefficient matrices and used 
during most of the node interaction computa t ions . 

The physical properties of the channel are 
calculated next in subroutine CHANAL. Here the north and 
South latitude boundaries, which were pre-set in INITGH by 
the user, are used to comipute the grid spacing along the x 
and y axis. Since this channel simulates a belt around the 
earth, the magnitudes of both EELTAX and EELTAY (meters) are 



proportional to the width of the channel divided hy the 
nur'cer of ^rid points in the y axis. 

The Cartesian coordinates for each node are ccrrputed 
'ty suhrcutine LCCaTI using the DILTAX and EILTAY calculated 
in CEANAL. If the option to use varying grid geometry is 
desired, subroutine TRANS transforms the grid geometry. 
TRANS also computes the corresponding new Cartesian 
coordinate values for each grid point and calculates the 
minimum IEITaX and TELTaY within the domain. When the 
geometry is changed to create a smaller ESLTAX or TELTAY, 
the two dimensional advective stability criteria is also 
changed. A new time step DT has to be computed using 
equation III-26. Since TRANS transforms the geometry, it 
also computes the new ET . 

Another transformation is required as discussed in 
Chapter II. The transformation from Cartesian coordinates to 
area coordinates is needed to perform the area integration 
of the inner products. Sub-routine AREA- computes these 
transformations exactly as outlined in Chapter II» Section 
C. Again cyclic continuity is very important and special 
care is needed to insure proper transformation. 

Ecllowing the area transformation is the computation 
of all the inner products that are required to solve the 
equations. Yhe advantage of using area coordinates is that 
the inner products (function of space coordinates only) are 
computed and stored once and used repeatedly without 
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rec2lcv.lation . Subroutine INN2H computes and stores these 
products using the formulas derived in Chapter II. 

The coefficient matrix, dimensioned NxN , where N is 
the totsl number of nodes, is a matrix of coefficients whose 
rows are the equations of the system to be solved. As 
discussed in the computational technique section, the 
members of this sparse matrix are the inner products for the 
left hand sides of the equations. Three coefficient matrices 
are used in the solution of the equations. The diagnostic 
equations (III-7 through III-15) use a coefficient matrix 
with the inner product <V . ,V, > which is constructed by 
subroutine Am, THU and stored in compacted form in vector 
GfN’N) by subroutine ASZm3L. However, when solving the 
prosnostic equations, these coefficient matrices have a DT 
''time step in seconds) term, so that these matrices are not 
assembled until the time integration begins. The vorticity 
and divergence equations (111-24,111-25) use the coefficient 
matrix H(N’N) with inner products 

solving the Poisson equations for the stream function and 
velocity potential, respectively. The continuity equation 
(III-23) uses a combination of inner products in its 
coefficient matrix J(NN) as follows 






$( At) 



2<Vj,Vi> 



IV-1 
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to solve the Helrrholtz equatioa. At the start of each tiire 

integration irodule, subroutine AMTRX2 is called to construct 

the two coefficient rratrices H and F. 

These banded and sparse rratrices are ccrrpacted into 

vectors to save storage during their assemblage by A3I^'3L. 

The vectors are dimensioned NN, as is NA^:E, the connectivity 

vector, and both use ISTART and NUM to index through therr. 

This compaction routine was used by Kelley (1976) and Older 

'1951) in their models, but was developed b 2 /’ Hinsman (1975). 

To illustrate matrix assemblage using an element by 

element technique, consider Figure 8. Note that this 

illustration is for element number 3, but all elements are 

treated in a similar iraner. The com.putational technique 

required that for each point (node) describing element 3, 

namely nodes 2, 3 and 14 stored in array FLMENT, the inner 

nro-uct <V.,V > between those points be distributed to their 
j i 

proper location in the coefficient matrix. 

Subroutine AMTRXl builds the lnner= product nodal- 
Interaction and stores it in matrix 3, dimensioned 3x3. 
Fig’.re S illustrates the 3 matrix for element 3, where the 
inner product <V . ,V, > is the multiplicand of the 
corresponding basis and test functions, respectively. 

The local dispensing of interactions is done in 
ASZ''-BL. Consider the second row of [3] in Figure 8. These 
are the interactions between node 3 of element 3 to the test 
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Figure 



AMTHX - builds [, 3 ] for ore element, passes pj and the element number 
associated with [ 3 ] to ASEKBL. This process continues till all element 
node interactions are assembled in the coefficient matrix. 



[Bl, = 



'V2 



^ V 

2^3 


"2"l4 


' V 

3 3 


"3"l^ 


14^3 


’'14^14 



ELMKBR 



ASEMBL - assembles the node interactions into the coefficient matrix 
[c] , which has the same structure as NAhll. The following diagram assembles 
inner product 5-r.to the coefficient matrix [c] , for element J. 
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14 


3 


15 1 1^* 
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144 
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pseudo-code 
•do 1-1 — 3 

11 - SI«EIiT(l) 



2 

3 

do j - START(li) 
jj . !TAME(j) 
do k - 1 -♦ 3 

kk - 2LM2hT(k) 



UST(ll) 



2 

3 

14 



1 




5 


6 




5 




1 


-»• -1 




16 




N 


21 


\ 


5 


26 






5 


31 






5 










1035 








1040 






5 






151 

152 

153 

15^ 

155 

156 



If ( kk - jj ) then 
C(j) - C(j) ♦ B(l.k) 
else continue 
Lend do 
-end do 
Lend do 



NAME 



10 

STAHT(3) -»11 
13 

last( 3) -.15 
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13 
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15 
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4 


3 


15 


16 
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155 


144 


1^5 



C(l) 



C(2) 



C(3) 



:(io) 



c(ii) 






10 

11 



cU3): ^ 



C(1U) 



c(i5) 



c(i6) 



C(20) 



15 

16 



20 



C(10h2) 



C(l0h3) 



C(lCij4) 



1043 

1044 



8. Assembling and storing the coefficient matrix for element 3. 



62 



functicns. ASZ^IEL locates nodes 3's connectivity list in 
using ISTART and NUM. In Figure S, this list is 
delineated by START(3) and LAST(3). Now ASEMBL steps through 
the connectivity list for three iterations. Turing each 
pass. ASI.'^EL is searching for one of the three node nuirhers 
for element 3. When a match is found with one of element 3's 
nodes (i.e. 2,3 or 14) and node 3's connectivity list (i.e. 
3.2,14,15 or 4) the proper position, to which this 
interaction is to he added has been found in the coefficient 
matrix. Since MAf^E and vector C, the compacted coefficient 
matrix, are dimensioned identically, the same pointer (i.e. 
J in Figure 3) is used to index through both arrays. This 
procedure is repeated for all elements in the domain to 
assemble the coefficient matrix of the equations. The 
pseudo-code for ASEMBL is shown in Figure 8 to facilitate 
stepping through this example. 

The domain and all finite element work vectors are 
initialized at this point. Subroutine ERb!S2T. is .called- -latter, 
to compute interpolation points for the harmonic analysis 
subrcutines . 

The last phase of the initialization process is the 
initialization of the dependent variables. The three input 
fields geopotential heights, stream function and velocity 
potential are computed in subroutine IC using the equation 
set III-30. However, the variables calculated from the 
diagnostic equations have to be computed using the input 
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fields, -hese varieties are used during each tirre step while 
solving the prognostic equations. 

The diagnostic equations are solved in subroutine 
BZ?Vn?., first during the initialization phase and later 
during the time integration phase. Each diagnostic equation 
calls its own module to compute the value of the forcing 
function and stores the computed values in the vector RES. 
These equations all use the same coefficient matrix when 
solving the diagnostic equations. Subroutine SOLVER is 
sufficiently genereal to solve each equation. SOLVER uses 
vector RES and coefficient matrix G to under-relax towards 
the solution. As mentioned previously, the coefficient 
matrix is strongly diagonally dominant so that three passes 
over the domain are sufficient. At the end of LSFVAR, output 
is generated depending on what print, plot, or harmonic 
analysis controls were preset. 

This completes the initialization pnase of the model 
and the program for the forecast phase will he described 
next . 



3 . Forecast Phase 

The forecast phase is accomplished in two steps. The 
first time step is made using two half steps by subroutine 
f^ATZNO. Here the prognostic coefficient matrices are 
constructed using half the LT value by calling AMTRX2. 
A^TRXE uses the same computational technique to construct 
the coefficient matrices as described for AMTRXl. 
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Zach of the prognostic eouations HI-23, III-24 and 
III-25 calls its cwn subroutine (CONTIC, VORTEQ and riVEC 
respectively' to compute all the terms on the right hand 
side, which are stored in the vector HHS. After computations 
for HHS are completed, subroutine RELAX solves the equations 
by cver-releicticn as described in the computational 
technique in Chapter 3. Once the solutions for the (n+1) 
time step are completed, DIRAR is called to update the 
variables from the diagnostic equations. Two oasses through 
YATZNO advances the solution fields one time step. 

The remainder of the forecast period is integrated 
using the leapfrog scheme. Subroutine LIAPFR performs this 
integration using the identical format as '^ATZNO, except 
that LT equals two DT . At preset times as specified in 
IN'ITGH, the different fields are saved for printing. This 
process continues until the final forecast time is reached. 

5. UTILITY yoruLis 

Once the equation formulation is completed, as in 

Chapter III, all the inner products and types of 

integrations are hnovn. Versatile modules can oe written to 

perform these computations. Consider a term of general form 

<^A . V . ,V . > where i is the node about which the term is 

evaluated and the j's are the nodes connected to node i, or 

the surrounding nodes. The inner product values <^V.,V. > are 

J ^ 

already computed and stored for all the nodes, during the 
initialization chase of the model. 
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-he reTdining coirputa ti on to complete tne evaluation of 
this term is the multiplication of the scalar coefficient of 
A at node J vith the corresponding inner product <V.,V.> for 
node This requires indexing through node i 's connectivity 
list stored in vector NAmz, and for each node in the list 
multiply and total the products. The cumulative sum of these 
multiplications is assigned as node i 's contribution for 
this term. Subroutine TIHm3 performs this exact computation. 
All that is passed to TERM3 is the scalar field A and the 
sign of the inner product, TZRM3 then computes the 
contribution for each node in the domain and accumulates it 
in the work vector ?.ES . 

Three other utility modules are> T2RM1, which computes 
the first scalar multiplication for triple inner products 
(i.e. ^Aj Vj Vjj ,V^ > ) . The product is again already 

computed and stored by subroutine INNER. TERMl computes 
Vfc .7 i "> and constructs a compacted vector similar to the 
coe-ficient matrices. This- reduces the ef-fort of- multiplying 
the second scalar to a TERf^3 computation. TZRf^2 computes 
node interaction of the following type <Aj 7ix^ * where both 
the basis and the test functions are derivatives. Lastly, 
subroutine TIRM4 computes node interaction for terms as <A.V. 

J JX 

,V^. , where only the basis function is a derivative. 

Vhen examining the right hand side of the equation sets 
III-23, III-24 and III-25, it is obvious this implementation 
is a subscripting nightmare; however, the use of the utility 
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rr'O'iules TIHMl, TIRM2 , TSRM3 and TERM4 siirplifies the 
irpleirentatlon to only determining what order to call the 
utility modules. Ixamination of subroutines CONTEQ, VCRTEQ 
and EIVIC, which compute the right hand sides for III-23, 
III-24 and III-25 respectively, illustrates this fact. No 
other subroutines or calculations were required. 
I’^pl emen ta ti on of these equations required minimal effort. 
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V. PHIMITIVZ MOLZL EXPZHIMINT 



The previous two chapters presented the detailed 
for’^ulctior and iirplen^entation of the vorticity-di vergence , 
shallow water equation rrodel. The results frorr this rrodel 
will he corrpared with the results from the primitive model 
in Chanter VII. To facilitate interpretation of the 
comparisons, a brief description of the primitive model 
follows. See Selley (1976) for a detailed discussion of the 
entire model. 

Also presented in this chapter is an experiment which 
demonstrates significant Improvement of the solution from 
the primitive model. Kelley's implementation used elements 
which were right triangles. Older (1S81) showed that 
eouilateral elements are far superior to triangular 
elements. - This experiment re implements the primitive model 
using equilateral elements and a comparison is made between 
the results of both implementations. 

A. .'"CTZL rZSCF.IPTION 

A form of the barctropic, shallow water, primitive 
equations developed by Phillips (1959) is used as the 
governing equation set for this model. In Cartesian 
coordinates the equation set is 
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height 
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velocity 


components at 


the same 



nodal points. This is an important consideration, iDecause 
the other models (the linearized model, see Chapter VI, and 
the vortici ty-divergeace model, see Chapter III) either 
stagger the dependent variables or have the property of a 
staggered formulation. Vhen com.parisons are made between the 
models, it is this lattice structure that is being compared. 

This form of the shallow water equations Includes 
gravity waves as a solution. Gravity waves have a maximum 
phase speed of about 300 meters/second. When the correct 
time step is calculated using equation III-26, a 
considerably smaller time step is obtained compared to the 
larger time step permitted in the vort ic Ity-divergence 
formulation. This is an important feature. If solutions from 
all models are equally as good, the best formulation would 
be determined using the computational time required to 
produce the desired forecast. 
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All Todels use the same domain structure. In fact, the 
dcrrain described in Chapter III was patterned after the 
domain implemented by Kelley. A,^:ain, this domain simulates a 
belt around the earth, with cyclic continuity which 
eliT^inates the east-west boundaries. Rigid boundary walls 
exist at the equator and at 30 degrees north latitude. The 
domain is composed of a 12x12 point mesh and subdivided into 
the right triangular elements illustrated in Figure 9. 
Notice that the grid points are not shifted as in the 
equilateral element implementation shown in Figure 6. 

The following boundary conditions are imposed: 

1) - no cross channel flow at the latitude 

boundaries . 

2) - a geostrophic balance at the channel walls 

imposed on the continuity equation V-1. 

This model has a simple second order diffusive term in 
the equations of motion V-2 and V-3. However, for the 

purpose -of evaluating these different-formulations, this 
option was not im.plemented during the comparison phase. 

Initial conditions consist of a single wave in the x 
direction and a half wave in the y direction. The initial 
fields for the three dependent variables are shown in Figure 
12. The maiimum zonal wind perturbation of 5.5 meters/second 
is superimposed on a mean zonal flow of 10 meters/second. 
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Figure 9 . Domain divided into right triangles. 
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Figure 10. Initial fields for the primitive model. Both the x and y 
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axes are multiplied oy 10 Km. 
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This Cursory description of the primitive model mentions 
only those sls-nificant features that will weigh heavily in 
the ccmparisions later. The Galerkin implementation of this 
model is similar to that presented in Chapters II and III 
and the system of equations is solved using a Gauss-Seidel 
iterative procedure. Further details concerning this 
pri'^itive model are given in Selley ( 1976). 



5. P.ZSULTS 

This experiment involves the shape of the elements. As 
mentioned previously, Kelley's implementation subdivides the 
domain into right triangular elements, as illustrated in 
Figure 9. Considerable small-scale noise was observed by 
Selley in the 43 hour forecast solution. 

The transition from right triangles to equilateral 
triangles changes the size of the domain. With right 
triangles, the and^y grid spacings are equal (300 KM). A 
12x12 grid matrix has a length and width of '3600 KM. With 
equilateral triangles, the^i and ^y grid spacings are no 
longer equal. Arbitrarily, the Ay grid spacing is held 
constant (300 KM) and a newAX grid spacing computed by 

4X = A.y/cos(30) V-4 
A 12x12 grid matrix with equilateral elements has a width of 
3600 KM and a length of 4045 KM. 

Figure 11 contains the 48 hour forecasts produced using 
both types of elements. The three dependent variables fields 
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Figure 11. ^4^ hour forecast comparison from the primitive model using 

both right triangles and equilateral triangles for element shapes. 
Arbitrary perturbation amplitude of 5*5 m/s. 
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are coppared for each. The small-scale noise that Selley 
observed is present. The three fields shew varying degrees 
c: distortions, which are especially noticeable along the 
boundaries . 

Cider (1981) found that the^“Oot trean square error (Hf^SE) 
was reduced 22 percent by using equilateral shaped elements. 
This improvement is apparent on viewing Figures 11b, d and 
f. The contours are smooth and the boundaries are 
noise free. Eelley showed excellent treatment of wave 
propagation by this primitive model. The lowest resolution 
grid '6x6) tested by Kelley was within four percent of the 
true phase velocity. Changing the element shape had no 
anparent effect on the phase velocities. 

Pecause the outcome cf this experiment was a 
significantly improved forecast solution, future comparisons 
with the primitive model will be made using equilateral 
elements . 
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Pigxire 12 . Different shaped basis fionctions . 
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A. ZCUATION P.ZTORMULATICN 

Zhe prirltive forr of the shallow water equations 
oresented In Chapter V (V-l,V-2 and V-3) is used to derive 
the linearized equations needed for this eiperirrent . The 
velocity equations V-2 and V-3 retrain unchanged and a linear 
basis function ^ ) is used to approiirrate the u and v 
var iabl es . 

The continuity equation V-l is linearized as follows: 




negligible 

where 'f is the average geopotential over the domain. A 
piecewise constant basis function (V . ) is used to 
approximate the geopotential. This linearization is 
reasonable in this case because the Ross’oy radius of 
deformation I /f* is much larger than ax [see Villiams and 
Zienkiswicz (1981)]. The. Galerkin methods is applied" -to- thi s 
linearized equation set using a piecewise linear test 
function for V-2 and V-3, and a piecewise constant test 
function for VI-1. 

The oiecewise constant basis function has the property 
of displacing the geopotential to the centroid location of 



the elements. 


which 


shoul d 


give 


the 


same effect 


as 


staggering the 


grid 


points. 


as 


seen 


in Figure 13. 


The 



density of geopotential data in the domain is now greater 
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Figiire 13 . Staggering of the geopotential 
about the velocity. The o symbol are the 
actual grid point locations and the ♦ symbol 
the centroid location of each element. 
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than the density of the velocity data. Instead of a 
gecpctential value for each node, there is one averaged 
value ever each element. 

The final forir of the Galerhin equation for VI-1, V-2 

and V-3 after perfortring the tire differencing is: 



^ J - v''<V VI-2 

- Vk^’rV-’'i- - ’'=-3 



VI-4 



This linearized set of equations VI-2, VI-3 and VI-4 is 
solved using a lauss Seidel iterative procedure. It is worth 
rentioning that the coefficient matrix <W. ,V,> in equation 

V 



V I -£: 



:as all non zero coefficients equal to one, since the 



integration of the inner product involves piecewise 

constant functions. 

This equation set is implemented rather than the 
equation set V-1, V-2 and V-3 using all of the existing 
pri*"iti7e model code. The major modification involved the 
way that the geopotential was referenced. With an average 
geopotential over each element instead of a value at each 
node fer a 12x12 mesh, there are 2S8 geopotential 
versus the lc5 velocity (u,v) points. 



oc i n t s 



B. HBSULTS 

The results from this linearized model are compared to 
those from the primitive model. The initial field for this 
experiment, Figure 14 has a maximium perturbation zonal wind 
that is one fifth of the value used in Chapter V, Figure 10. 
The mean zonal wind remains 10 meters/second. The 48 hour 
forecast solutions for each field are compared in Figure 15. 

This alternative formulation shows some promise, 
although there are some minor perturbations in these 
contours compared to those in the primitive solution. No 
explanation is offered for this small-scale noise, although 
possibly the other formulations presented by Williams and 
Zienkiewicz (1981) would improve the solutions. 
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Figure 14. Initial fields for both the primitive model and the 
linearized model. 
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Figure 15 . ^8 hour comparison Between the primitive model and the 

Linearized model. (aPA = 1.1 m/s) 
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VII. 7C?.TIClTT-riVZRGZ.\CE MCZ'IL ZXPIRI^!ZNT 

•-.s iz the previous two chapters, a coirparisor between 
two models will be presented. The results from the 
VO r t i ci tv-d i ve rgenoe model are compared to the results from^ 
the primitive model. To fully exploit the differences 
between the models performances and indicate the strengths 
and weaknesses of each formulation, three domain geometries 
and three initial conditions are used. All solutions are at 
4S hour, except for one case which was extended to 96 hours. 

Zrom these two finite element formulations some 
additional insight is obtained concerning the execution time 
rec.uired as the grid resolution changes. Lastly, a brief 
discussion on the sensitivity of the computational technique 
is ^'iven. 

A. TZST TOMAINS AN'L INITIAL CONIITIONS 

The three domain geometries used in the model evaluation 
are illustrated in Figure 16. All domains consist of a 12x12 
element mesh with equilateral shaped elements (156 grin 
points' and cyclic continuity is imposed on the east and 
vest boundaries. The domain has dimensions of 4345 KM along 
the X axis and 3533 KM along the y axis. 

The regular domain (Figure 16a) has a uniform 
distribution of grid points, with a minimum grid spaci-ng 
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Figure l6 . Test domain geometries . 
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alorg the x axis of 337 XM . This is the Tost oongenial 
domain ar.d orcduces the test results. 



The s’^octh dorair. ("igure 13’o) has a smoothly varying 
distriouticn of grid points. This technique allows a smooth 
variation of resolution. It was developed by Older (19S1). 
who showed how it significantly reduced noise at all 
-requencies compared with other variable grid domains. 

The degree of resolution variation is accomplished by 
selecting an appropriate value for the ratio of maximum 
stret^'h to minimum shrink along both axes. The experiments 
presented in this chapter use a 2.c ratio for both 
directions. This produces a grid point concentration in the 
right center of Figure 16b and coarse resolution at the top 
and bottom left of the domain. The minimum grid spacing 
along the x axis is 199 SM in the dense grid point area. 

The third domain was the least hospitable geometry for 
both models. The abrupt domain (Figure Icc) has a dense grid 
point concentration on the left and coarse spacin? on the 
right of the domain. The minimum grid spacing along the x 
axis is 15S KM in the area on the left. The grid spacings 
are uniform exceut for the abrupt change along the center. 

Although these three domains are simple in structure, 
they are adequate to test both models' performance. To 
further enhance some di f f erentiating characteristics between 



the two formulations, three initial conditions are used, 



wo 



have previously been described in Chapters V and VI. Their 



S6 



rain ii stinsiuishin^ feature is the artitrary perturta ti cn 
arulituce ra^nitude. The nearly linear case has an APA 
= 1.1 .t/s ccrpared tc 5.5 r/s for the more nonlinear case. 
All initial conditions have a 12 m/s rean flew component. 

A'ith the nearly linear case both rodels behave well. The 
introduction of the more nonlinear initial disturbance 
illustrated the boundary and computational technique 
sensitivity. The third initial condition is the nearly 
linear initial field with the wave length equal to half the 
domain length. This has the effect of producing two waves 
with the domain . 



2. T-ST CASI COb'PAHISCNS 

The comparisons between models are divided according to 
the domain geometry. The primitive model has three fields? 
geopotential, u and v. The vcrticity — divergence model has 
seven fields: geopotential, stream function, velocity 

potential, u, v, vorticity and divergence. The geopo ten tial , 
u and V will be the only fields used for the comparison. 

1 . Regular Case 

The regular domain geometry (see Figure 16a) is used 
•■’or this first comparison. The initial conditions have an 
A?A =1.1 T/s, as shown in the contour plots in Figure 17. 

The 4:3 hour solutions for both models are presented 
in Figure IS. As anticinated, all of the solutions have 
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Figure 1?. Initial fields for both the primitive and vorticity — 
divergence models using the regular domain and perturbation amplitude 



of 1.1 m/s . 
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Figure 18 . ^8 hour forecast comparison between the primitive model 
and the vorticity-divergence model using the regular domain. (aPA = 
1.1 m/s ) . 
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srocthiv 'ontoured fields with no noticeble noise and their 
phase velocities are comparable. 



■''ith the nearly identical solutions, the 
d i 5 1 i nsui shi r.^ feature between the models is the 
compute ticnai time. The computational time is determined by 
the size of the time step each model is allowed to use as 
indicated on ecv.ation III -2c. For this domain the grid 
spacing is 337 The maximum phase velocity possible is 
differert for each model. The primitive model allows gravity 
waves so c - ZZls m/s whereas the vortici ty-di mergence model 
filters out these high frequency waves. A c = 10 m/s is used 
for this formulation. 

Table 1 



*"0061 


ox 


c 

( m/s ) 


<:!^t 

( sec ) 


of steps 
(iterations) 


CPU units 
( sec) 


? ” 


337 


300 


455 


375 


145.5 


V-n 


337 


10 


13,753 


13 


33.0 



Comoarisons of the computational times for a 4B hour 
forecast between the primitive model and the 
vorticity-divergence model. 



Table 1 co^’pares the results of the computational 
times for both models. The vorticity divergence model 



produced as accurate results over the regular domai.n \:Sing 
22% of the computational time needed by the primitive model. 
In fact from geostrcphic reasoning the vort icity-divergence 
model should produce better forecasts when small scale 
features are oresent. 
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Sroo th Case 



The sjooth dorain ^eorretrv (Tig’jre idh) is used in 
this second comparison. The i.nitiai fields shown in Ji^rure 
Ir have a disturbance amplitude of 1.1 r/s. 

The 4S hour solution comparison is presented in 
rleure 22. All fields a,§ain have smooth contoured plots. 
Close inspection of the two u fields, ri,?urs5 20c and d, 
shows that the nrimitive u field has small kinks along some 
contours and a weak tilt near the central boundary ncdes, 
althou^rh this may be a function of the plotting routine, 
h'otice the good symmetry for the vorticity-divergence u 
field . 

As in the regular case, this smooth experiment 
produced two acceptable solutions. Again, the computational 
tim® is the differentiating criteria. 









Table 2 






'"odel 


AX 


c 


A t 


a of steps 


CPU units 






(m/s ) 


(sec ) 


(iterations) 


( sec ) 




ISG 


300 


271 


c38 


304. S 


V-D 


199 


10 


8124 


21 


43 . £ 


cm par 


i so n of 


the corrpu 


t a t i 0 n a 1 t i 


mes for a 45 


hour forecast 


etwee 
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^ 1 ."m i t i V s 


model and 


the V 0 r t i c i 


t y-dive rgence 


od el 


using the smooth 


doma i n . 






able 


2 comuares the 


computational times for 


bo th model s . 


he VO 


r t i c i t y 


-divergence has a 33. 


8% saving of 


CPU time. 
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Figure 19. Initial fields for both the primitive and vorticity- 
divergence models using the smooth domain and perturbation amplitude 
of 1.1 m/s . 
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Figure 20. ^8 hour forecast comparison between the primitive model 
and the vorticity-divergence model using the smooth domain. (aPA = 
1.1 m/s) 
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effort to contrast the co.rputat icnal accu 
;his experin'ent is repeated using the 
:ase. Figure 21 shews the initial fi 
rr/s . This larger initial disturhanc 
reflected in the greater geopotential ar'plitude 
ragnitudes of the contour lines. 



f the models. 


this 


on 1 i near initial 


case . 


ith A? A = 5.5 


m/s . 



racy 

Tore 

elds 

sis 

and 



Figure 22 shows the 42 hour solution cotr.pari son . As 
in the rore linear case presented above, all of the plots 
have smooth contours with no noticeable noise. The 
vort ici ty-di vergence geopotential field, Figure 22b, has the 
ridge extending farther north, and flattening of the 
southernmost contour, than does the primitive geopotential 
field. Figure 22a. The mean geopotential heights for both 
models have also increased. The primitive mean geopotential 
is now 492S2 gom and the vor tied ty-d ivergence geopotential 
is 49620 gpm . 

These two discrepancies indicate that the boundaries 
are not handled accurately. Turing the implementation of the 
vortici ty-di vergence model, treatment of the boundaries was 
the most troublesome phase. The vor tici ty-di vergence 
formulation is a complex equation set and time limitations 
restricted further investigation of more sophisticated 
boundary conditions. 

This same initial condition is now extended to a 96 
hour forecast, which is shown in Figure 23. The mean 
VC r t i c i ty-d i vergen ce ^eopotential increased to 49S20 gpm. 
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Figure 21 . Initial fields for both the primitive and the vorticity- 
divergence models using the smooth domain and perturbation amplitude 
of 5-5 ni/s. 
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Figure 22. ^8 hour forecast 

and the vorticity-divergence 



comparison between the 
model using the smooth 



primitive model 
domain. (aPA = 5*5 



m/s) 
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Figure 23- 96 hour forecast comparison between the primitive model 
and the vorticity-divergence model using the smooth domain. (aPA = 
5.5 m/s) 
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vhereas the rrean priiritive geopotential retrained constant. 
All fields have sn^ooth contours. Close inspection of both u 
fields, figure ?3c and d, shows a slight skewing of the 
contours along the central channel grid points. The 
vorticitY divergence u field has the ir.ost pronounced 
deviations. A hypothesis for this skewing, explained further 
below, is that it is caused by the corputa t ional technique 
employed. The relaxation schemes are extremely sensitive and 
fine tuning of the relaxation coefficient would have 
required more time than was available. 



Table 3 



f^od el 


A X 

(T.y.) 


c 

(m/s) 


a t 

(sec) 


m of steps 
( iterations ) 


CPU uni 
( sec) 


?I 


199 


320 


271 


1276 


554.2 


V-D 


195 


12 


8124 


42 


91.2 



Ccrnarison of the computational times for a 9€ hour forecast 
between the primitive model and the vortici ty-divergence 
■^odel using the smooth domain. 



Table 3 shows the comparison between both models'for 
the 96 hour forecast. A savings of 55 percent is realized 
with the vorticity divergence model. 

The above experiment points out the two areas where 
the vorticity — divergence model is presently weak, the 
increase of the geopotential and the sensitivity of the 
relaxation coefficients. 3oth these weaknesses can be 
i'^proved and are not a result of the formulation but of the 
i mo 1 emen t a t i 0 n . At oresent, their influence is not detected 



o T - 0-3 t in 

- e"o n ' t ra * e 



extrernely long forecasts, such as in the 96 hour 
rther experirrents ray still be required if they 
significant differences in the solutions c: the 



tvc red els. 

The last experirent on the srocth dorain uses the 
rcre linear initial condition APA =1.1 r/s, but its wave 
length is divided by two, so that two shorter waves are 
propagated through the channel. The initial conditions are 
shown in Pigure 24. Eecreasing the wave length has the sare 
effect as decreasing the density of grid points. In this 
case six grid points are used to describe the wave structure 
versus the 12 used in the previous cases. 

The 46 hour corparisons are shown in Figure 25. As 
in all previous cases, a computational time saving of S4% is 
^^ained with the vcrtici ty-divergence model. With fewer grid 
points describing the wave structure, more small-scale noise 
is introduced into the solution. Comparing the primitive 
geopotential field, Figure 25a, to the initial geopotential, 
Jisure 24a, shows a dampening of the wave amplitude, whereas 
the VO rtici ty-divergence geopotential. Figure 25b, 
correlated well with the initial geopotential. 

The high frequency noise is evident on both u 
fields. Figures 25c and d. The primitive model u field is 
poorly defined along the boundary and becomes irregular over 
the interior grid points. 
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Figure 2^. Initial fields for both the primitive and the vorticity- 
divergence models using the smooth domain. There are two waves embeded 
in the flow and A?A = 1.1 m/s . 
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Figure 25. ^ hour forecast comparison hetween the primitive model 

and the vorticity-divergence model using the smooth domain. Two waves 
are emceded in the flow and APA = 1.1 m/s . 
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The smooth dorraln allows a variatle resolution of 
points and produces excellent results. As the wave 
ler.'th s’ets shorter, or the forecast length gets larger, 
scTe small-scale noise is apparent, especially with the 
primitive formulation. 

3 . Ahruot Case 

The atrupt case comparison uses the abrupt domain 
geometry ^see Figure 16c). This grid point configuration is 
used to •'‘urther illustrate the effects cf spatial 
resolution. The previous case using the smooth domain and 
half wave length introduced noise into the sclution, hut the 
spatial resolution changed slowly and gradually. 

Consider the transition necessary in an operational 
model, where the luxury of' haviug a-iong' smooth-'transition- 
into the region of high resolution may not he possible. The 
atrupt domain is an example of the results obtained when 
spatial resolution is decreased rapidly. 

The initial fields are shown in Figure 26. The more 
linear case, A?A = 1.1 m/s, is used to eliminate effects due 
to the initial field, so that only the effects due to the 
eri'’ geometry are seen. 

The 4£ hour comparisons are shown in Figure 27. 3oth 
solutions are affected by this geom.etry, hut the primitive 
solutions are totally disorganized and unacceptable. 

Table 4 shows the comparison of computational times 
for both models for a 48 hour forecast. An 35S savings in 
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Figure 26. Initial fields for both the primitive and vorticity- 
divergence models using the abrupt domain and perturbation amplitude 
of 1.1 m/s. Note that the element shapes are not equilateral triangles. 
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Figure 27 • ^ hour forecast comparison between the primitive model 
and the vorticity-divergence model using the abrupt domain. (APA = 
1 . 1 m/s ) . 
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'--"'U tire is obtained using the vorticity-divergence model, 
yet only is there a computational savings with the 
vort i ei ty-di ve rgen-'e model. hut the solution is 
significantly better than the primitive model. 

Table 4 



Mo del 


J X 


c 


At 


^ of steps 


CPU units 






(m/s ) 


( sec ) 


(iterations) 


{ sec) 


p7 


16? 


330 


228 


756 


368 . 9 


* i- 


168 


10 


6858 


25 


55.7 



Comparison of the computaional times for a 48 hour forecast 
between the primitive model and the vcrticity-livergence 
model using the abrupt domain. 

C. COMPUTATIONAL SINSITIVITT 

This section will offer an e.rplanation for the skewing 
o*' the contours In the 96 huur forecast suiut ion* over -the- 
smooth domain 'Pigure 23). As mientloned ureviously, there 
was not enox’.gh time available for fine tuning the 
overrelaxation coefficient, which is used while solving the 
system of eouations. The overrel axati on coefficient may be 
very sensitive and small changes can, on occasion, radically 
charge the rate of convergence. The optimal value of the 
c verrelaxa ti on coefficient depends on the specific ferm of 
the coefficients of the equation and the error distribution. 

The ecuation set to be solved consists of three 
eouations and each equation required its own relaxation 
coefficient. >‘hen solving the equations over the regular 
domain, the entire system is well behaved and an optimum 
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relaxation coefficient is easily detertrined. Eowever, as the 
dor-aln ?ecretry changes, the system cf eouaticns dc not 

converge as ranidly and the relaxation coefficient requires 
other fine tuning. 

The 96 hour forecast uses the strooth domain. The 
mid-la ti tvdinal grid points are compacted, creates a denser 
celt in the middle of the channel. The coefficients 

originally computed using the regular domain need 

adjustments to properly solve the equations. 

To illustrate the significance for fine tuning the 

relaxation coefficient, consider the series cf plots in 
Figure 23. The vorticity divergence equation set can he 

simplified hy assuming the flov is non-divergent , so that 

only the vorticity equation needs to he solved. Figure 23a 
is the 48 hour vorticity field using this equation over the 
regular domain with an o verrel a.xa tion coefficient of 1.3. 
The field is well defined with smooth contours. 

Figure 28h is similar to the case in Figure 28a except 
that the smooth domain is used. Notice the V-shaped kinh in 
the oattern with a steeper slope in the upper half. This 
i'’creased bias in the upper half is caused hy relaxing the 
field in the same direction during each pass over the 
domain, '/.'hen the direction is reversed after each pass, the 
exaggerated bias in the upper half disappears, as is seen in 
Figure 2Sc . However, the V-shaped kink is not eliminated. 
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Figure 28. Computational sensitivity using the ^ hour vorticity fields. 
Fig. 23 a) the ^ hour forecast using the regular domain, 28 B) the same 
forecast using the smooth domain relaxing in one direction, 28 C) same as 
28 3) except alternate the direction of relaocation after each pass over 
the domain, and 28 D) sajne as 28 C) except the relaxation coefficient was 
fine tuned form 1.3 to 1.297* 
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Varying the relaxation coefficient from 1.2 to 1.2S7 
crocv.ces the nrnch iriprcved solution in Figure 28d . If the 
relaxation coefficients for the other two equations could 
also he fine tuned, it appears that improved solutions would 
result. There are also other relaxation techniques available 
that have potential for improving the solution while also 
converging at a faster rate. Some of these techniques will 
be tested in the futuer using this model. 
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v::i. CONCLUSIONS 



This research investigated different finite element 
f orrrulations for the shallow water equations. The 
twr — dirensional dorrain was a channel which sixulated a telt 
around the earth. Analytic Initial conditions were used to 
siv'plifY the ccxparisons. Two forxulations were exax.ined; 
one using different shaped basis functions and the ether 
using a different forx of the equation. Zach was compared to 
the orixitive form of the shallow water equations that was 
developed by Selley (1975). 

The use of equilateral shaped elements which was 
suggested by :r. w.J.F. Cullen significantly improved the 
solutions compared to lelley's model, which originally used 
right triangles as basis functions. .'*^ost of the other 
studies in this thesis used the equilateral triangles. 

'Williams and Zienkiewicz (1981) suggested the use of 
piecewise linear basis functions for the velocity field and 
piecewise constant functions for the height field. This 
formulation was tested with a linearized continuity 
equation. The results were poorer than those obtained with 
lelley 's model . 

yost of the effort in this thesis was devoted to 
ixplementating and testing a vcrticity-divergence model 
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similar to the ones developed by Stanifortb and i^itchell 
(19'’'?) and Cullen and Hall (1979). Several tests vere 
presented which corrpare this forrrulation with Kelley's 
Tcdel. It was found that this mcdel executes approxi.Ta tely 
one order of magnitude faster than does the primitive 
foT'Tulstion used by Kelley. Secondly, as the spatial 
resolution between grid points decreases, this formulation 
oroduces a solution that is far superior to the primitive 
form. A disadvantage is its computational sensitivity, which 
requires fine tuning in solving the elliptic equations for 
certain geometries. It also requires 25 percent more 
computer storage, due to the more complex equation set and 
the additional variables that are treated. 

Implementation of finite element miodels is not easy. 
However, there is a lot of generality and redundancy 
i-bedded in the computations. Versatile modules were written 
which si gni f i cantly reduced the effort in implementing the 
vor tici ty-d ivergence model. 

rurther research is suggested using this finite element 
formulation. It has accurate phase pr onagation , is able to 
handle variable grid geometry, reduce the small-scale noise 
and decrease the model's execution time. Specifically more 
advanced methods of solving the elliptic equations should be 
investigated. Tinally, the formulation should be tested with 
small-scale forcing, where its advantages should be most 
evident . 
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